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Abstract 

The Gaussian beam superposition method is an asymptotic method for computing high 
frequency wave fields in smoothly varying inhomogeneous media. In this paper we study the 
accuracy of the Gaussian beam superposition method and derive error estimates related to 
the discretization of the superposition integral and the Taylor expansion of the phase and 
amplitude off the center of the beam. We show that in the case of odd order beams, the error 
is smaller than a simple analysis would indicate because of error cancellation effects between 
the beams. Since the cancellation happens only when odd order beams are used, there is 
no remarkable gain in using even order beams. Moreover, applying the error estimate to the 
problem with constant speed of propagation, we show that in this case the local beam width 
is not a good indicator of accuracy, and there is no direct relation between the error and the 
beam width. We present numerical examples to verify the error estimates. 

Keywords: wave propagation, high frequency, asymptotic approximation, Gaussian beam superposition, 
accuracy, error estimates 

1 Introduction 

Simulation of wave propagation is expensive when the frequency of the waves is high. In this case, 
a large number of grid points are needed to resolve the wave oscillations, and the computational 
cost to maintain constant accuracy grows algebraically with the frequency. At sufficiently high 
frequencies, therefore, direct simulations are no longer feasible. 

Instead one can use high frequency asymptotic models for wave propagation. The most popular 
one is geometrical optics, which is obtained when the frequency tends to infinity. The unknowns 
in geometrical optics are the phase and amplitude which are independent of the frequency and 
vary on a much coarser scale than the full wave solution. They can therefore be computed at 
a computational cost independent of the frequency. However, a main drawback of geometrical 
optics is that the model breaks down at caustics, where geometrical optics rays intersect and the 
predicted amplitude is unbounded. 

Gaussian beams approximation is another high frequency asymptotic model which is valid also 
at caustics. It was introduced by Popov [T], based on an earlier work of Babic and Pankratova [5]. 
A Gaussian beam is an approximate high frequency solution to the linear wave equation which 
is concentrated close to a standard ray of geometrical optics, called the central ray of the beam. 
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Although the phase function is real-valued along the central ray, Gaussian beams accept complex- 
valued phase functions off their central ray. The imaginary part of the phase is chosen such that the 
solution decays exponentially away from the central ray, maintaining a Gaussian-shaped profile. 
The main advantage of this construction is that it gives the correct solution also at caustics. It 
has recently been proved to be beneficial in seismic imaging, (6117]. 

Numerical methods based on Gaussian beams use the superposition principle. Individual beams 
are computed via ray tracing like equations, where quantities such as the curvature and width of 
beams are calculated from ordinary differential equations (ODEs) along the central rays, and 
the contribution of the beams concentrated close to their central rays are determined by Taylor 
expansion. The full wave field is then obtained by a superposition integral over all beams. This 
integral is replaced by a discrete summation of beams in practical computations. See for example 
[31 HI [5] ini [Z] • Numerical techniques based on both Lagrangian and Eulerian formulations of the 
problem have been devised [51 IHl UHl El HI] • For a rigorous mathematical analysis of Gaussian 
beams we refer to [T5] . 

In this paper we derive error estimates for the beam superposition method. We study the 
discretization error, caused by replacing the superposition integral by the summation of beams, 
and the error related to Taylor expansion of the phase and amplitude off the center of the beam. 
Some error estimates for this method have been derived earlier, jl4| , I15 j . We aim to give a more 
complete picture of the error by also including the error due to the spreading of the beams, which 
is related to the Taylor expansions. This error is recognized as important in e.g. [14 ] . It turns out 
that, in the case of using odd order beams, the error is smaller than a simple analysis would indicate 
because of error cancellation effects between the beams. Since the cancellation happens only when 
odd order beams are used, there is no remarkable gain in using even order beams. Moreover, 
we show that in the case of constant coefficient equations, i.e. when the speed of propagation is 
constant, the local beam width is not a good indicator of accuracy, and there is no direct relation 
between the error and the beams' width. However, this may not be true in the case of varying 
speed of propagation, where the beam width can be an important factor in the Taylor expansion 
error. For other recent results on error estimates see [l6l ITT]. 

In Section 2, we review the construction of Gaussian beams and the Gaussian beam superpo- 
sition method. The accuracy of Gaussian beam superposition is studied in Section 3, where the 
main result is formulated together with numerical examples verifying the obtained error estimates. 
In Section 4, the proof of the main theorem is given in detail. Finally, in Section 5, we compute 
the errors analytically in the case of constant coefficient equations and give some remarks on how 
to select the Gaussian beam parameters. 

2 Gaussian beam superposition method 

Gaussian beams are obtained when the linear wave equation is solved with oscillatory initial or 
boundary data with an amplitude in the shape of a Gaussian bell. A Gaussian beam is an asymp- 
totic solution concentrated on its central ray in the domain. By the superposition principle for 
linear equations, such solutions can be added to find the full wave field. The initial/boundary 
data for beams are obtained such that the wave data at the source is well approximated. In this 
section, we consider the Helmholtz equation and review the construction of Gaussian beams and 
their superposition. 
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2.1 Construction of Gaussian beams 

Consider the Helmholtz equation 

^2 

Au(x) + —-ru(x) ^0, xencR^, 

C[X)'' 

where a; ^ 1 and c{x) are the frequency and speed of propagation, respectively. Boundary con- 
ditions are given on dfl, which we assume is divided in two parts: one where ingoing waves are 
specified, and one with outgoing radiation condition, typically at infinity. We call the first, ingoing, 
part of dH, the source curve. We substitute the WKBJ ansatz 

oo 
fe=0 

into the Helmholtz equation. Here, the phase function and the amplitude functions Ak are 
assumed to be smooth and independent of lo. Equating coefficients of powers of lo to zero gives us 
the eikonal equation and the transport equation for the phase and the first amplitude term in the 
frequency domain, 

\V4)\^l/c{x), 2VAo-V0 + .4oA0 = O. 
For the remaining amplitude terms, we get additional transport equations 

2 VAfe+i • V0 + Ak+i A0 + AAk = 0. 

When u! is large, only the first terms in the WKBJ expansion are significant. We henceforth denote 
the high frequency approximation taking K terms in ([T]) by uqo{x)^ 

UGo{x) = A{x)e'"^^^\ A{x) ^ Ak{x){iu;)-^ (2) 

fc=0 

This approximation is usually called the geometrical optics approximation, in particular when 
K ^ 1. It introduces an error of the order 0{uj^^). 

The Gaussian beam approximation has the same form as the geometrical optics approximation, 

ugb{x) ^ A{x)e'^^^^\ (3) 

where the phase (j) and amplitude terms Ak satisfy the same PDEs. There are, however, two 
important differences. First, while in geometrical optics is globally defined for all rays, for 
Gaussian beams it is constructed based on one specific ray (the beam's central ray). Second, in 
geometrical optics, <p is real- valued, but in the Gaussian beam construction it is real- valued only on 
the central ray of the beam. Away from the central ray, it is complex- valued with positive imaginary 
part. The solution will then be exponentially decreasing away from the central ray, maintaining 
its Gaussian shape. Note that since cj) is complex valued, it actually satisfy the complex eikonal 
equation, \18\ I19j . Unfortunately, the question of existence and uniqueness of the complex eikonal 
equation is to a certain extent still open. In particular what precise boundary conditions are 
well-posed for the above setting is not known. 

As in geometrical optics, the Gaussian beam approximation breaks down when (/^{x) becomes 
non-smooth. This is typical for solutions to both the standard and the complex eikonal equation. 
It happens in general some distance away from the central beam. On the other hand, away from 
the beam the solution rapidly goes to zero and the precise value of the phase is not important. 
One usually deals with this problem by multiplying the amplitude with a smooth cut-off function 
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that is one close to the central ray, and zero for some fixed distance away from it. In practice, ([3]) 
is thus replaced by 

where ^p{x) is smooth and compactly supported around the central ray. 

For a beam starting at point Xq with direction Pq, the corresponding central ray satisfies the 
ray tracing ODEs 

dx , , , dp Vcix) , ^ , ^ Pn 

at at c(x) \pq\c[xo) 

with t being the real- valued travel time along the ray. If we set p — {cos6, sin 6')^/c(a;) and 
X = {x, y)^ we can reduce (|4]) to 

dx dy / N . d9 , ^ . „ , ^ ^ 

— =c a; cost^, — =:c(a; smt^, — = c^. a; sm — c„ a; cos 6'. (5) 

dt ^ ' dt ^ ' dt ^ ' ' ^ ' 

The complex- valued Ak and (j) close to the central ray are then approximated by Taylor expansions 
around the ray, 

Ak{x) « Ak{x*) + {x- X*) ■ S/Ak{x*) + ^{x - x*yD^Ak{x*) {x - x*) + ■ ■ ■ , (6) 

0(a;) w (P{x*) + {x- X*) ■ W<j>{x*) + ^{x - x*y D^(t>{x*) {x-x*) + --- , (7) 

where x* — x{t) for some t. The Taylor coefficients (j){x{t)), V(j){x{t)), Ak{x{t))^ etc. on the 
central ray can be computed. The lowest ones are real on the beam and given directly, 

^{x{t))^^{xo)+t, Wc^{x{t))^p{t). 

The higher order ones can be obtained by solving ODEs similar to Q, and may have a complex 
part on the beam. The most common approximation by far is to take K — I in ([21), so that 
A{x) = Aq{x), and to approximate A{x) to zeroth order and 0(a;) to second order. In this case 
we have, [S], 

A{x(t)) = A{x,) [^-^ ^) , D^c^{x{t)) ^ HNH^, (8) 

with 

H ^ ( ''"^r ^^^f) , N^( n , h) ^ i/^Vc, (9) 

y-cos0 sm^y y—ci/c — C2/c^/ \C2 J 

and the complex-valued scalar functions P and Q satisfy the dynamic ray tracing ODEs 

^=c\x)P, Q{0) = Qo, (10) 
dP Cx^sm^9-2cxySm9cos9 + CyyCos^e 

The quantities P and Q determine the leading order wavefront curvature and the beam width. For 
example, if y = a; — a;* is orthogonal to the beam at x*, then by (jSj and ^ 



\uG^{x* +y)\ 



showing that the effective beam width is proportional to [uj'^{P/Q)/2\ It can be proved that 
if Qo 7^ and 9(Po/Qo) > 0, then Q{t) ^ and '^{P{t)/Q{t)) > along the central ray for all 
t > 0, [1]. Therefore, by a proper choice of initial data Qo and Pq, each beam will be regular (with 
finite amplitude at caustics) and concentrate along the central ray. A common choice is Qq > 
and Pq = i. 
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Figure 1: The sum of several Gaussian functions is almost constant. A plane wave can therefore 
be decomposed approximately to a sum of parallel Gaussian beams. 



2.2 Beam superposition 

Let the source curve be given by xo{s) in parameterized by s. We introduce the notation 
A(x,s), (j)(x,s) and ^p{x,s) for the amplitude, phase and cut-off of a beam with initial position 
Xo{s). In the Gaussian beam superposition method, the boundary condition on Xq[s) for the 
wave field is asymptotically expanded into Gaussian beams, [5l. Individual Gaussian beams are 
computed by solving the ODEs (jlj and (llOllip . The contributions of the beams concentrated close 
to their central rays are determined by the approximations ([6l [7]) entered in ([3]) . The wave field is 
then obtained by the superposition integral over the beams. 

In practical computations, this integral is replaced by a discrete sum of individual beams, the 
trapezoidal rule approximation, 

u'^{x)=uj^''^h'^ip{x,Sj)A{x,Sj)e"^''''^''^'^\ (13) 

where h is the initial spacing of the beams. 

The initial conditions for the Taylor coefficient ODEs are chosen such that itf well approximates 
the exact ingoing boundary data. This can be done in different ways. In particular the initial 
width of the beams can be varied to give different approximations. As an example, we consider 
a plane wave in the y-direction as boundary condition on the x-axis, xq{s) = (s,0). This will be 
approximated by a sum of beams starting in the same direction, [6] . The approximation is based 
on the relationship 

1 = f e-^--'^"'''^'ds = V ^ + 0(e-('"'/'')'), = jh, (14) 

with 770 representing the initial beam widths, see Figure[TJ Identifying (fT4|) with ([T2l[T3| . assuming 
= 1 we see that 

1 {^-sf 
A(x^ 0, s) = : — , 0(a;,O, s)=« 5 — . 

To properly choose the initial data, one must take the parameters 770 and h such that rjQ > h 
by ([13]). Then the wave field ((T^ will produce an accurate plane wave on Xq{s) = (s,0). The 
condition rjo > h can be related to the initial data {Po,Qq) of the dynamic ray tracing ODEs 
(|10llip . Since beams go in the j/-direction 6 — tt/2 and we have H — I and 4>xx — P/Q by ^ and 
®. Thus, 

Pq _ 2i 
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Chosing Pq = i we therefore get 



h<r]o= 



With this relation between h and Qq we get an accurate approximation. In particular we need a 
spacing h of order 0{l/^yuj) for a fixed Qo- Note also that for computational efficiency, h should 
not be taken much smaller. These restrictions were derived for a plane wave but similar scalings 
will be necessary also for more general boundary data. 

In what follows, in order to simplify the calculations, we assume that all beams, originating 
from Xo{s), shoot out orthogonally. We denote hy X{t,s) the location of the center ray originating 
in xo{s) after time t. We further assume that (j){xo{s),s) = 0. 

We make one observation that will be used in the analysis bcilow. It is well-known that 
Xt II Va;!/), Xt ■ Vx</> = 1 and Xg -L X^ under the assumptions made above. Therefore, since 
<l){X{t,s),s) = t, 

= -^HXit, s), s) = X, • + c^,= M^{t, s), s) (15) 
as 

Hence <ps = everywhere on the central rays. 



3 Accuracy of Gaussian beams summation 

In this section we study the accuracy of summation of Gaussian beams. One can distinguish six 
different types of errors in the approximation: 

1. High frequency approximation. 

2. Error in initial data. 

3. Discretization error. 

4. Taylor expansion error. 

5. Cut-off error. 

6. Error in numerical integrators for solving Taylor coefficient ODEs. 

The first error depends on the number of terms used in the WKBJ approximation, i.e. the difference 
u{x) — Us{x). For example, for standard beams it is of the order 0{l/ui) since one amplitude term 
is used, meaning that each beam is a solution to the Helmholtz equation up to order 0{l/uj). The 
second error represents how well the exact boundary data is approximated by a superposition of 
Gaussian beams. The third error is caused by replacing the superposition integral by a discrete 
summation of beams, i.e. Us{x) — uf{x). The fourth error is due to the fact that A and (f) 
are not computed globally, and only their derivatives on the central beams are computed. One 
therefore needs to approximate their values around the central beams by Taylor expansions. The 
fifth error is caused by multiplying the solution by a smooth cut-off function in order to account 
for possible irregularities away from the central rays. Finally, the last error is the numerical error 
in solving the ODEs for computing Taylor coefficients. For example the global error in a fourth 
order Runge-Kutta method is O(At^), with At being the time-step. 

Here, we will only concentrate on the Discretization and Taylor expansion errors. 
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3.1 Motivation and preliminaries 



Let the source be a curve Xq{s) and assume that we look for the solution along a Ime x — {x,y*). 
We simplify the notation by setting Ak{x,s) = Ak{x,s), 4>{x,s) — s), Lp{x,s) = ip{x,s) and 
write 

M,(a;)=a;i/2 f ^p(x, s) A{x, s) e"^'^^''-'^ ds, 



.D 



(x) = uj^/^hj^fix, Sj)Aix, sj)e*"*("- 



Sj = jh. 



We now let X{s) denote the location of the center beam on the line (a:, y*) when the initial data 
is given at Xo(s). Hence, X{s) is implicitly defined by 

Xitis),s) = {X{s),y*), 
for some function t{s). Figure [2] shows the setting for Xo{s) — (s,0), as an example. 

y 



X = {x,y*) 


Xis) 


Xo{s) ^ (s,0) 





Figure 2: A schematic representation of the initial source and a beam central ray. 

We will now explain the approximation of A{x, s) and (f>{x, s) for a, {q + l)-th order Gaussian 
beam, complying with the standard notation that the basic choice g = is a first order beam. We 
then take K — lq/2\ + 1 terms in the WKBJ expansion ([2]) and observe that the high frequency 
approximation error will be of the order 



0(a;-«'/2) where q* = 2([9/2j + 1) = 



(7 + 2, q even. 



q + 1, q odd. 

For the term Ak{x, s) we make a Taylor expansion up to order q — 2k around X{s), 



(16) 



Ak{x, s) « Ak.q-2k{x, s) := AkiX{s),s) + ■■■ + 



^J^^^^^^dr'^A,{X{s),s), (17) 



with < fc < [q/2\ . We also set 



L9/2J 



Aq{x,s):^ ^ Ak^q-2k{x,s){iuj) ''^A{x,s). 



(18) 



fc=0 
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Furthermore, we approximate 4){x, s) up to level g + 2, 



The approximate Gaussian beam solution is then given by 

u,ix) ^^^'^ j v{x,s)Aq{x,s)e"^^''^''^'Us, 

The reason why the phase is approximated to two orders higher than the amplitude is to balance the 
Taylor expansion errors; the phase error is multiplied by the frequency uj, which is proportional to 
one over the beam width squared (cf. pO| below). Note that for q> 2, one needs to take K > I in 
([2]), i.e. to include more terms in the WKBJ expansion in order to also balance the high frequency 
approximation error and the Taylor expansion error, cf. Remark [2] below and the discussion in 

m- 

Our motivation for considering the Taylor expansion error comes from the following observation. 
We define the width of the Gaussian beam passing through {x,y*) as 

T]{x) := . 

y^UJ '^(j)xx{x,X-^{x)) 

Because of the term ^^^i'^-^i'^)) "^-xx/a ^^le solution will be close to zero for \x — X{s)\ > rj{x). A 
simple error analysis would therefore give the following result. Using ^ we have 

Us-Us = {A- Aq)e"^^-' + ^e*'^<#-,(e»"('^-'^.) - i) 

= {iio)-''{Ak - ifc,g-2fc)e'"'^' + ^e*'^^9(e^'^W-'?<!) _ l) 

A;=0 
L9/2J 

fc=0 

Hence, since 77 — 0{(jJ^^/'^) we would have 

k/2j 

Us-u^^Y 0{io-^'i+^y^) + ^e^"*'.0(c^-(9+3)/2+i) ^ 0(c^-(9+i)/2). (20) 

In particular, for first order beams with g = 0, the convergence rate in w would be half order, 
i.e. proportional to l/^/uj. This is also what is observed numerically for a single Gaussian beam. 
However, we now consider two numerical examples using superposition of first order beams to 
verify the convergence rate for this case. Since it is difficult to obtain the exact Gaussian beam 
superposition solution Us we here instead compare Us with first order {K = 1) geometrical optics 
UQo, which is close enough to Us to verify or refute the sharpness of ()20p : the high-frequency error 
in both Ms and first order geometrical optics is of the order 0(l/w). Hence, the predicted error of 
first order beams, 0(1/^/0;), would dominate if it was sharp since 

\y-s - "Gol < \Us ~ Us\ + \Us - u| + |u - ugoI < C{l/^/uj+ I / Uj) < C/^/uj. 
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In the first example, a plane wave generated on the line y = propagates orthogonally into the 
computational domain with a variable speed of propagation. Figure [3^ shows the central rays of 
Gaussian beams, and Figure shows the absolute value of the Gaussian beams and geometrical 
optics solutions along the line y = 0.6, shown in bold in Figure[3^. Figure|3l3 shows the logarithmic 
scale of the maximum error between the Gaussian beams solution and the geometrical optics 
solution. As can be seen, the convergence rate of the error is surprisingly proportional to u!~^, 
which is half order better than what we expected. 

In the second example, a plane wave generated on the line x = propagates with an angle of 
45° into the computational domain with a variable speed of propagation. The convergence rate of 
the error, shown in Figure [3f, is again proportional to lo^^. 

We will therefore study the Taylor expansion and discretization errors more carefully to describe 
why it is smaller than what we expected. 

3.2 Main result 

For our results we make the following precise assumptions 

(Al) Smoothness of all coefficients. We assume Ak{x,s) e C^'^'^~''^~^'^(IR^), the space of functions 
with p + q + 2 — 2k continuous and bounded derivatives. Similarly (j){x, s) G (jp+i+'t ^nd 
Xit,s) e CP+i, withp > 2. 

(A2) Algebraic growth of phase off center beam. For pi,p2 < P + 9 + 4, we have for some p, 

dP^dP^^{x,s)<Cil + \x~Xis)n. 

In particular, all derivatives are bounded on the center beam, x = X{s). 

(A3) No caustics. The derivative X'(s) is bounded away from zero, < cq < X'{s) < ci < oo. 

(A4) Non-degeneracy of each beam. The imaginary part of (jjxx on the beam is strictly positive 
and bounded, 

< Co < '^(l)xx{X{s), s) < ci < C50. (21) 

Moreover, the frequency is non- vanishing, lj > C2. This means that the approximate beams 
will have a fast decay off the central beam for high frequencies, and also that the beam width 
never vanishes or becomes infinite. The last point is an important feature of Gaussian beams, 
related to the fact that Gaussian beams can approximate the exact field at caustics. 

(A5) Cut-off of fixed size. We use (p{x,s) — (p{X{s) — x) with tp e C°° such that (p{x) = 1 for 
|a;| < a/2 and (p{x) = for |a:;| > a. The size of a will be chosen "small enough" depending 
on (f> but independent of uj. 

The error that we want to estimate is given by 

E{x) = Us{x) - uf{x) = Us{x) - Us{x) + Uis{x) ~ uf{x) ~. E'^ + E^, 

where E"^ — Us{x) — Us{x) and E^ — Us{x) — u^{x) represent the Taylor expansion error and the 
discretization error, respectively. 
Then we can show 

Theorem 1 (Main Theorem) For the (g + l)-th order Gaussian beams, we have 

\us{x)-uf{x)\<\E^\ + \E^\, 
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(e) (f) 

Figure 3: Left and right top figures show the central rays of Gaussian beams by an initial plane 
wave on x- and y-axis, respectively. Middle figures show the absolute value of the Gaussian beams 
and geometrical optics solutions along the lines y = 0.6 and y = 2. Bottom figures show the 
logarithmic scale of the maximum error between the Gaussian beams solution and the geometrical 
optics solution. The convergence rate of the maximum error is ui~^. 



where 




\E'-\<Cu;-^, q* = r:7 ZT ^ (22) 

q odd, 
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and 

|£-|<C(-A.)'^ (23) 

The constants depend on p, the initial data, Pq and Qq, for the beams but does not depend on x, 
uj or h. For the Taylor expansion error we have 



i.e., the leading order term of the error E"^ in uj is C*{x) ui^^'^ r]'^ ^ ui with C*{x) given 

by and 

Remark 1 // we take h < ri{x), the discretization error E^ is typically smaller than the Taylor 
expansion error iS"^ because of the "spectral" accuracy in i2S\) . For first order beams (with q = 
), the observed convergence rate is therefore first order in uj, which is the same as geometrical 
optics. However, h should not be chosen too small for computational complexity reasons. It is also 
important to note that to balance the error with the error in initial data, h should also relate to the 
initial beam width rjQ. 

Remark 2 The estimate shows that the Taylor expansion error indeed balances the high 

frequency approximation error I116\) . Moreover, it suggests that there is no remarkable gain in 
using even order beams (with an odd q); neither the high frequency nor the Taylor expansion error 
get a better convergence rate in l/io with these beams. However, one should note that this is only 
true in the case of the superposition of beams, where error in adjacent beams cancel. If we only 
have one beam, this does not hold and the simple error estimate in i2U\) is sharp. In this case the 
Taylor expansion error dominates the high frequency approximation error for even order beams. 



4 Proof of main result 

Before going on to the proof of Theorem [1] we prove the following utihty lemma concerning 
estimates for the composition of two functions. 

Lemma 1 Suppose gs{z) belongs to C^(R) for each value of the parameter S. If 

|gf <C^fc(l + |^r), l<k<p, (24) 

where Ck and q > are constants independent of z and 6, then there are functions h„i^k G C^^''"(M) 
and constants Cm,k independent of z and d, such that 



-e 



^ hr^^z), max \h^-\{z)\ < C„,,,(l + {zl"'). (25) 



dz'^ ^-^ " 0<n<p 

m— 1 



Proof: We show ((25|) by induction. For = 1 we have /ii.i = g'si^) ^ ^ and the statement 
clearly holds. Suppose ([25|) is true for 1 < fc < p. Then 



jfc+i k k+i 

__easiz) ^ ^ssiz) ^/^;„_,(z)+g^(z)/i„ ,(z) = e^^W ^/i™.;c+i(^). 

m— 1 m— 1 

Thus 

iKi,k^ m = l, 

hm,k+i{z) = I h'^^j. + g'ghm-i,k, I <m<k, 

yg'^hrn-i^k, m = fc + l. 
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Using the induction hypothesis, we immediately get that /im,A;+i(^) G ^ ^(^)- Moreover, 

n 

max |/.L".UiWI<,<^ax + max X'^n\hlil,Jz)U"+'-^\z)\ 

n<p—k~l ' 0<n<p—k—l ' 0<n<p— /c— 1 ^ — ' ' 



0<n<p-fe-l "'.'-T^-^ - 0<ri<p-fc-l 0<n<p-fe-l 

J=0 



The first term is bounded by Ci^i (1 + by assumption, and for tlie second term we can estimate 



wfiicfi proves (|25|) . □ 

We can now start witli tlie main proof. We consider each error separately. 

4.1 Taylor expansion error 

The Taylor expansion error is given by 



(x) - Us{x) = J <^(x(s) - x) s) e*'^'^^^''*) - s) e''^'^''^^''*)) ds (26) 

7 n 



k=0 

In this subsection we will start by studying the general integral approximation 



= J (p(x(s) - a;) s) e*"*^"'") - i,, {x, s) e*"^«("'^)) ds, (27) 

where, with a slight abuse of notation, A and will represent one of A^. and Ak q_2k respectively 
in the sum above. We can therefore also assume that Qa l£ q and that q — qa is even. 

Let us denote X^^{x) by m{x) and then, since X'(s) is bounded away from zero we can use 
the change of variables 

z= — — ^ s = m(x^^{x)z). (28) 

^(x) 

We obtain 

= ^1/2 77 y (p(77z) (a(x, m(a; + ryz)) e''^'!^(^^"'(^+''^))- 
(x, m(a: + ryz)) e*'^'^^(^^™(^+''^))) m'(a; + t^z) dz. 

Now, letting 

Da{x,s) := A(a;, s) - Aq^(x,s), D^{x,s) := (/)(a;,s) - if)q{x,s). 
and recalling that supp(^ C [—a, a], we can write the integral as 

E'^ = ^1/2 77 / 93 (Da + A(e^"^'^ - 1)) e^^^-'m'dz. (29) 

We will now approximate the terms in the integral (|29p by their Taylor expansion. Let us use 
the shorthand 

(— 1)P ~ 1 dP 

apix) = —dPA{x,m{x)), bp{x) = — — A{x , m{x + z)) 

p\ p\ dzP 



z=0 
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and 



^p(^) = —J-9x<i>ix,m{x)), fp{x) = - — (t){x,m{x + z)) 



pi 



p\ dzP 



(30) 



2=0 



We note that, in this notation 



q+2 



Let 



Ag^{x,m{x + z)) = '^dj{x + z) z\ 4>q{x,m{x + z)) = '^pj{x + z) z^ . 
j=o j=0 



ai{x) = dg^+i{x), a2{x) = dq^+2ix) + a'^^+i(a;), 



bi{x) 



ci{x) = 5R 



{x,m{x)) ' 
r2{x) 



b2{x) = i 



C2ix) 



.Pq+4{x) +p'g+3{x) 



Q(l)xx{x,m{x)) 

.hix) - (Tps{x) 



,{x,m{x))' " '<s(f)xx{x,m{x))' 

where cr = 1 for g = and cr = for > 0. We then approximate 

DA{x,m{x + 'nz)) w r)'^-+'^DA{x,z) := (r?^;)««+^ai(x) + (r?^;)««+^a2(a;), 
giu,D^(a.,m(x+r,z)) _ ^ ^ r]i+^B{x, z) := r]''^Hi{x)z'i^^ + r?''+2(62(a;)0«+* + ah\{x)z'^''+^ /2), 
^iui4>,{x,m{x+r]z)) ^ c{x,z) =: e'"''^(^'"*(^''+'^^'=i(^'~^^/^(l + C2{x)r]Z^). 

The residual terms are denoted 

Da{x, m{x + r]z)) - rj^^+^DAix, z) =: rj^^+^^RAix, z), 
^iu,D^(x..m{x+nz)) _ i_r,'^+iB{x,z) =: 7]'^+'^ Reix, z), 



ciu4>q {x,m(x+r]z)) 



C{x,z) =: TfRc{x,z). 



Then we have 



Lemma 2 Lei the residual terms Ra, Rb and Rc he defined as above. Under assumptions (Al) 

and (A2), for small enough a, 

\RA\<C\z\'i-+\ \RB\<Ce'''\ \Rc\<Ce-'"^\ V|0| < a/r?, 
where the constant C is independent of x, lo and z. 

Proof: Wc note that Uq^ {x + z) are the first Qa coefficients in the Taylor expansion of A{x + z — 
x', m{x + z)) around x' = 0. Therefore, by Taylor's theorem and assumption (Al) 

\Da{x, m{x + z)) - ^«»+^<+i(a; + z) - 0«»+2<+2(a; + ^)| < C\z\''-+^. 

Expanding the second and third terms around z = gives the bound for Ra- 

We now estimate coD^ in two different ways. By Taylor's theorem as above, for some (, with 

< rjz, 

\LuD^{x,mix + vz))\ < \di'>+^Ui^,mix + Tjz))\ ^^f^ < Cc^lry^l'+'Kl + |r?^n, 
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where we used the growth condition (A2) for to bound the error term. Then, for \z\ < a/rj, and 
small enough a, 



implying 



(31) 



Moreover, the same steps as for Da together with (A2) gives 

\D^{x, mix + z)) - z'^+%+s{x) - z'+Hp,+i{x) + P,+3i^))\ < ^^^^+'(1 + {zf). 
and since ui — l/?7^50xx, when \z\ < a/r], 

\iujD^{x,Tn{x + rjz)) - 7]'^+^ z'i+%{x) - r]'>+^ z'i+%{x)\ < CY+^\z\'^+^. 
Finally, for q > 0, clearly \ujD^\^ < Ctj'^''+^\z\^i+^ < Cr;9+3|z|29+6 and for q = we get 



\{tiuD^)^ - r,^blz^\ 



\Dl - {n^z^p3f\ _ _ r^^z^p^WD^ + r;3z3p3| 



< Ct^Iz 



Thus, 



\R^b\ < C7729|z|39+9e^'/« + ^ (1 _ a)C|z|29+6 + aC|z|29+7 < C'e/'^^. 



To show the third inequality, we note that since (f)s{x,m{x)) = by (llSp . we have fi{x) = 0. 
Therefore by Taylor's theorem and assumption (A2), for q' > 2, 



{x, m{x + z)) — (j){x, m{x)) — z^fp(x) 

p=2 



<qzr'+i(i + |zn 



(32) 



Let v{x, z) = 4>q{x, ra{x + z)) — 0(x, m{x)) — z^f2(x). Then, by ((3T|) and (|32]) . 

\v{x,z)\ = \(l){x,m{x + z)) - 4>{x,m{x)) - D^(x,m{x + z)) - z'^f2{x)\ < C|z|^(l + \zf). 
Moreover, 



\v{x, z) - z\h{x) + ap3(x))| < C\z\\l + \z\P). 



As above, if \z\ < a/r], 



yu:v{x^r,.) _ ^ _ r]Z^C2{x)\ < \iUJv{x, Tjz) - T]Z^C2{x)\ + ^\ujv\'^ 6^'^"^ 

< Cij\'nz\*{l + \t]z\P) + i \CUJ\7]Z\^{1 + |,yz|P)|%'^"l''"l'(l + l''^l'') 

< Ct]^\z\\1 + an + Crj^\zf{l + aP)2eC^^"(i+"n < Crj^e/-"/\ 
for small enough a. It remains to note that, since (j)s{x,m{x)) ~ '^(j)x{x,m{x)) = 0, 



"^h = irn'(x)^50s5 = ^m'{x)'^ 



d 

dx 
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which shows that iujrf'z'^f2 — iz^ci — Therefore, 



and the estimate for \Rc\ follows. □ 

We Taylor expand the remaining quantities in ([29]) and use the assumption (A5) to get 

Lp(iqz) w 1, 

A{x, m(x + ?7z)) « A{x^ z) :— A{x, ni{x)) + 77261 (x), 
m! {x + rjz) « m(a;, z) :— m'{x) + rjzm" {x). 

It is easy to see that the residual terms for these approximations can all be bounded by Crf'z^. 
Since these residual terms as well as Ra and Rb above all grow slower than exp(z^/4), we can 
replace the terms in the integral in ([29]) by their approximations and control the error by 0(77''"+'^), 



ABCrhdz 



<Cuj^/^r^r)^-+\ (33) 



Moreover, since the C(x, z) is exponentially small in 77 for \z\ > a/i] this estimate holds also when 
taking the integral over all of M. We can now compute the leading error terms in 77. The first one, 
ei = en + 77612, is 



DACrhdz 

= e*""*(^'™(^)) J z'^^+^aim' + T]z'^-+\a2m + aun" + aiC2m' z^)e"^''^'''-'^-'^ '^dz + 0(77^) 
= : e''^'^("''"("»(eii + 'yeis) + 0(77"), 

where 



and 



en = dq^+iaim', = a2m'dq^+2 + aim"dq^+2 + aiC2m' dq^+4^, 

d,ix)= 1^.(1 Peven, 

i [0, podd. 



(34) 



(35) 



with Np being a constant. We note that dp{x) = when p is odd and it is bounded in x when p 
is even. The second term is 62 = £21 + ^&22, 

ABCmdz === J [A{x, m{x)) + 77261 (x)] [6i(x)z«+3 + 77 (b2ix)z''+^ + ^blix)z^'^+^^ 

^ g^<^0(:r,m(x))Wci(:r)-2V2(l + 77C2 (a;)^^ ) (77l' (x) + 7]Zm"{x))dz 

= e*'^^(^'™(-)) J ^Abiz^+^n' + 7y(6i6iz«+V + A (b2Z^+^ + ^blz^''+''^ m' 

+ Abiz''+^C2rn' + Abiz'i+'^m"^ e"'"^^''^''^ '"^dz + 0{jf) 
=: e*"*("'"("»(e2i + 7/622) + 0(77^), 
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with 



= Abidq+am', 622 = bibidq+4,m' + A ^62C?g+4 + -^bld2q+ej m' + Abidq+GC2m' + Abidq+^m" . 

(36) 

To find an expression for the leading order error term we now have to consider four cases depending 
on Qa- First, if ga < g then the second term in (I33p is of the same order or smaller than the right 
hand side and we can disregard 62- Second, if Qa is even then en — 621 — since dp ^ for p odd, 
and wc gain an additional order in rj. Upon also noting that q — qa is even, we can therefore write 



where 



<? = 



+ 2, qa even, 
+ 1, qa odd. 



C(x) 



iuj<^{x^m{x)) 



'en + 621, (?a odd and q = (ja, 

612 + 622, (ja even and g = Qa, 

611, qa odd and qa < q, 

,612, qa even and qa < q. 



(37) 



Moreover, C{x) is independent of uj and h and can be bounded by a constant independent of x. 

We now go back to the full Taylor expansion error in and use the results that were 
obtained above for (j27l) . Clearly, all parameters and functions will depend on the term number fc, 
and we indicate this with a subscripted k. Since g — 2fc is even if and only if q is even, we have 
qk = go - 2fc = q* - 2fc with q* defined in Then, 



L9/2J 

E^- ^(*c^)-'=Cfc(x)c.i/2,7?-+i 

k=0 



The result therefore follows with 



L9/2J 



Y,i^^r''(.El-Ckix)u;''-n 



/2 „<!k+l 



k=0 



(38) 



L9/2J L9/2J 

A:=0 fc=0 



L?/2J _ L9/2J 

^*(^) = E («^)"'^fc(a^)^7*"'* = E (*^^^ 

k=Q k=a 



2\-k 



Ck{x). 



(39) 



Since turf' = 0(1), the leading order term of the error E'^ in is 



-9*/2 



4.2 Discretization error 

The discretization error is given by 

E^ = Usix)-uf{x)=uj^^^ f ^(X(s) -x)ig(x,s) e*""^'("'^'ds-wi/2/iE /(■?■)' 



with 

f{j)=cp{X{s,)-x)Aq{x,s,)e^-^^^^'^^\ 
for a fixed a;. The Poisson summation formula gives 
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where 



f{k)= / f{s)e-^''''''ds^ / ip{X{sh)-x)Ag{x,sh)e''^^-''^^'''''h~^'''''''ds 



Therefore 



fe#0 

Using the change of variables ([^5]) we obtain 



/» 



(40) 



We will now show that the integrand functions in (|40p are smooth, with bounded derivatives. 
Then the non-stationary phase lemma can be used to bound f{k) since the phase derivative m'ix) 
never vanishes. 

We need 

Lemma 3 Under assumptions (Al), (A2) and (A4-), for < £ < p and \z\ < a/rj with small 
enough a, 



d^ 

dz^ 
d 



Aq{x, m{x + ryz)) 

ioj(pq[x,7n{x+rjz)) 



dz^' 



(41) 
(42) 



The constants C and C are independent of of £, x, lo and z. 



Proof: For the first inequality we can consider the individual terms in the sum psp separately. 
They will each be of the form 

q-2k 

Ak^q-2k{x, m{x + rjz)) = ^ dj{x + rjz) (rjzY . 

3=0 



where 



(-1)^' ■ 
aj{x) = —^diAkix,m{x)). 



By assumption (Al), dj £ are bounded, for < £ < p, uniformly in a; and, after noting that 

\r]z\ < a and that \ri\ is bounded by a constant because of assumption (A4), the result (|4T|) follows. 
For £ — the second inequality is obtained by writing 

e^'^^^i^M^+n^-)) = cix, z) + jfRc, \C{x, z)\ = \l + 7]z^C2{x)\ e-'-^''^. 

Now since \rjz\ < a and 02(2;) grows algebraically by assumptions (A2), (A3), (A4), and since 77 is 
bounded by a constant, we have by Lemma [2l 



\e'^M=^M^+vz))\ < \C{x,z)\+r]^\Rc\ < Ce-^'/^ 



(43) 
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Now consider 1 < £ < p. Wc write 



q+2 



4>q{x,m{x + r]z)) ^^pj{x + r]z){r]zy , 
with Pj{x) defined in ([30)) . Then, since pf, +pi = by ([T5|). we have 



— (f)q{x, m{x + ?7z)) = ffpiZ + (rj-'^^ Z-' Pj{x + r/z) + jrj-' z-' ^Pj{x + rjz)^ 



9+2 



and therefore 



^ rn(a; + rjz)) 



q+2 



where 



7j + (j + 1 < j < g + 1, 79+2 := Pg+2- 

Since the phase derivatives are evaluated on a center beam, 7^ € are bounded, for < £ < p, 
uniformly in a; by assumption (A2) and we therefore have 



^ [iuj(j)q{x,m{x ^--qz)) 



<Q(l + |z|''+2), l<i<p. 



Thus, by Lcmma[l]with — iuj(j)q{x,m{x + rjz)) and 6 = 7], using (j25p and (j43p . the inequahty 
(H^ fohows for 1 < £ < p. This completes the proof. □ 

The remaining terms in (j40p . i.e. ip{r]z) and r7i'(x + rjz), are all assumed to be smooth with 
derivatives of order up to p bounded uniformly in x by the assumptions (Al) and (A5). Since the 
growth in ([4T|) is offset by the rapid decay in p2)) . the above Lemma shows that all z-derivatives 
of the integrand, 

g{x,z) ■.= ^Aq e'^^- m', 

up to orderp belongs to Li and \\d^g{x, •)|| < Ck for Q <k <p. The constants Ck are independent 
of X and uj. We can then use the following version of the non-stationary phase lemma. 



Lemma 4 Suppose jj){z) e Cp+^(M) with € C^(K) and ^p'{z) > cq > 0. Moreover, let 

e < S < 1 and suppose g{z) € W^'^ . Then 



<C\\g\\w.^{^-^ , 



(44) 



where C depends on tp^x) and p, but not on g{z), 6 and e. 

For the proof we refer to [2^. It is an easy adaptation of the proof of theorem 7.7.1. 
Taking ij} as 27rm(x + •), 5 as r] and e as h/k we can apply this to 



I/>)I = I 



rj 



krj J 



Consequently, 



h \ri 



k-p < C 



k=l 



h \ j] 
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Thus since by the assumptions (Al) p > 2, 




Together with this shows the theorem. 

5 Constant coefficient equations 

It is often claimed that the beam width is important in the accuracy of Gaussian beams, because 
for wide beams the Taylor expansion error should be large. See for example [H [3] . We therefore 
in this section consider the constant coefficient Helmholtz equation, with the speed of propagation 
c{x) = 1, for which exact Gaussian beam solutions and the Taylor expansion error \E'^\ can be 
computed. We investigate the importance of the beam width on Taylor error in this particular 
case. Our conclusion is that the local beam width is not a good indicator of accuracy, and there 
is no direct relation between the error and the beams' width. We show the main steps of the 
derivation and the final expression for C*{x) and the leading relative error terms below. For more 
details we refer to [T2] . 

We consider first order beams where q = 0. These only contain one amplitude term Ao{x) which 
we for simplicity call A{x) here. The source curve will be denoted xo{s) = (s, yo{s)) and we assume 
all beams originating from Xq shoot out orthogonally. Therefore 9o{s) — ^ + tan^^(j/Q(s)). In the 
constant coefficient case the central ray is a straight line. With a;(0) = xo{s) = s, y{0) — yo{s) 
and 9{0) = 9o{s), we get from ^ at y ^ y* , 

e{t{s)) = ^ + Un-\y',is)), (45) 
x{t{s)) = X(s) = s - y'ois) [y* - yo{s)), (46) 
t{s) = {{X{s)~sY + {y*-y,{s)f)"\ (47) 

Here we will only compute the error at a; = (0,y*). For this point, let s* := to(0) = X^^{Q). To 
simplify the calculations, and without loss of generality, we assume yo{s*) — ?/o(s*) 0. Therefore, 
by (|45p - ((Tf|) . the central ray starting at Xo{s*) will lie on the y-axis, and we have s* = X{s*) — 
and t{s*) = y*. See Figure HI 



y 



He 

y = y 


iX{s*),y*) 







Figure 4: A schematic representation of the initial source and central beam rays which are straight 
lines. 
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Assuming the initial phase on Xo{s) to be zero, 0(a;o) = 0, we also get 

<j>iX{s),s)=tis), 0(O,TO(O))=y*. (48) 

To obtain ODEs for higher order Taylor coefBcients, we introduce the orthogonal ray-centered 
coordinates t,n, where n is the axis perpendicular to the ray at point t with the origin on the ray. 
In this coordinate system, (j)(t,n = 0) and A{t,n = 0) correspond to (/)(X(s),s) and A{X{s),s) 
in the Cartesian coordinate, respectively. The eikonal equation and transport equation in the 
ray-centered coordinates read 

(49) 

2VA • v./) -f AA(/) = 0, v./) = (0t (^„)^. (50) 

Set (t) di(j){t, n = 0) and A'-^^ (t) :^ d^^A{t, n = 0), with j = 0, 1, 2, . . . . We first note that 
by®, 

0W(i)=t, 5t0(i,n = O) = 1, a,Xt,n = 0) = 0, j = 2,3,.... 
Moreover, by (|49p and ([SO)) and taking several of their derivatives with respect to t and n, 

(j)W(t) = 0, ata„0(t, n - 0) - 0, dtdf,(l){t,n^O) ^-(j)^^^\t), 
dtdl(j){t, n = 0) = 0, d^d.a(t){t, n = 0) = 0, dfdnC^it, n = 0) = 0, 
52a2<^(t,n - 0) = 20(2)3 (i), ^^^^^^ „ ^ 0) = -^A^'Ht) <t>^^\t), 
d^A{t,n = 0) = ^A(o)(t)(^(2)2^^^^ atc)„^(t,7i = 0) = 0. 



Now, let 



and 



</)(i,n) « i + y0(^^(t) + -^^''(t) + ^^''\tl (51) 

A{t, n) « (i) + nA(i) (t) + y v4(2) (t). (52) 
Putting ([^ and ([5^ into (H^ and (|5(I1) . we obtain the following ODEs for the Taylor coefRcients, 

^0(2)+ 0(2)^=0, 



+ 10(2)^(0) ^ 0, 



d 



+ -0(2)^(1) + 10(3)^(0) ^ 0, 

dt 2^ 2^ 



+ ^0(2)^(2) + 20(3)^(1) + 1^(4)^1(0) + ^^(2)3^(0) 

diti 2 2 2 

We then solve these ODEs with ^("^(0) — 1 and zero initial conditions for the rest of Taylor 
coefficients. At our observation point x = (0, y*), we have that =9^, since the n-axis is parallel 
to the a:-axis. We can therefore easily transform the solutions in the ray-centered coordinates t, n 
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back to the coordinate system x, s. In the end we note that all terms with odd x-derivatives are 
zero. Hence, we obtain that ai(0) = 61 (0) = and ei2 + 622 in ([51)1 and ([55)1 simplifies to 

ei2(0) + 622(0) = m'(0) 02(0) d2(0) + m'(0) A{0, 0) 62(0) ^4(0). 

After some additional algebraic manipulations and assuming that Pq — i, SQo = Oi ^Qo > 0, we 
get 

"2(0) = « A(n ^ ■ *\7/2 ' 

, ... . {Ql + y*') + 4(Qo + »?/*)^m-(0) 

m^^TJ^^ ' ^''^ 

,^^o^^r + (Ql + yymoir)^ (55) 



mO)= .^ .../2 ^ '^W^f^^^is^) ■ (56) 



and 

A(n r\\ — ,,,,,, _ , 

Moreover, by (|45ll47p . 

^J{y*) = y'^iO), m'(0) = {X-')'{0) = (1 - 2/*2/^'(0))-^ . (57) 

Therefore, knowing yo{s) and by (|53ll57p and (|35|) . we can calculate 

ei2(0) +622(0) = e-''^<t>i0M0))c*[0) ^ e-'"^*C*(0). 

Note that C*(0) only depends on Qo, y* and yoiO). 
We now consider the following two canonical cases: 

(1) y'm = o, 

(2) y'^io) = -1. 

The first case corresponds to a line t/o = 0. The second case corresponds to a circle yo{s) = 
— 1 + Vl — or a parabola 2/0(5) = —s'^/2. Note that with an initial curve with positive second 
derivative, the rays will intersect and form a caustic, and then our theory does not hold. 
For the first case, we obtain the simple expression 

where no is a constant complex number. For the second case, the expression is much more com- 
plicated. In the small and large Qo-limit we have 

^*/n^ -iLjy' "^i + '^2y* + nsy*'^ 2 , rnfr,3\ r^*fn\ -i'->y' "4 ^-5/2 , rn!r>-'^/'^\ 

<-c 0) = 6 ,4/11 * Qq + 0{Qq), = e — ===Qo + 0{Qa , 

y*^y/l + y* Vl + y 

(59) 

where rij, with j = 1, ... ,4, are constant complex numbers. The amplitude of the geometrical 
optics solution is proportional to |1 — 

y* yo(0)r^^^ and by the relative error will be 
= \E'^\ |1 - y* 2/i'(0)|i/2 + 0{u-^'^) = u^^'Wi^) \C*m \l-y* y'm"^ + 0{u-^'^). 
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4 6 



(a) 



(b) 




1 1.5 




(c) 



(d) 



Figure 5: Absolute value of relative error as a function of Qq (top) and of the width -q (bottom) 
in the case when the initial source is a line (left) and a circle (right). 



We therefore obtain the leading order term 



\E, 



rcll 



and for small and large Qq, 

-1 ni 



corresponding to ([58]) and (f59l) . respectively. 

Figure [5] shows the absolute values of the relative errors at y* = 3. As it can be seen from 
the formulas and figures, the relative error has a direct relation with Qq, tending to zero both for 
small and large Qq. A reduced error with large Qq has also been noticed in [14] for the oscillatory 
part of the error (or the discretization error). However, there is no clear connection between the 
error and the beam width; the same width can correspond to different errors. 
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Figure 6: The beam width as a function of Qo at y* = 3. 



In many approximations, the optimal Qq, corresponding to the minimum beam width at a 
receiver point is chosen for computations, see [1] for instance. Figure [5] shows the beam width as a 
function of Qo- In our case the minimum width is attained at Qo = y* ■ With Qo = U* ^-nd y* ^ 1 
we get 

with N and N' being constant numbers. When using this Qq, we do not obtain the minimum 
error as was seen above. However, importantly, the relative error decreases as the distance from 
the source increases. 

We conclude that in the case analysed here large and very small Qo will improve the Taylor 
expansion error. From Figure |6] we see that this corresponds to having wide beams, not narrow 
beams. One should keep in mind, however, that this is not the whole story. The approximation 
of the initial data where the source curve is not flat and/or the amplitude is not constant will in 
general deteriorate when wider beams are used. Hence, this restricts the beam widths that can be 
used. Wider beams also mean that the wave field will be more expensive to evaluate since beams 
contribute more globally to the solution. Moreover, our result is strictly for constant coefficients. 
In the presence of a varying speed of propagation where the properties may change dramatically 
as we get farther from the central rays, the Taylor expansion error could be large for wide beams. 
In addition, when the rays can bend, it may not be possible to have very wide beams, since as 
was noted before, the Gaussian beam approximation may break down when the phase becomes 
non-smooth, and this happens at some distance away from the central ray (outside the regularity 
region). In the general case, finding the optimal Qo for a given observation point is an open 
problem. 
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